Intro

Main analysis: This scripts reads spatially referenced linear trends over time in biomass, metabolic index and temperature (“05_calculate_trends.Rmd”), and fits spatial sdmTMB models to those

# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(RCurl)
library(devtools)
library(sdmTMB)
library(rMR)
library(patchwork)
library(raster)
library(VoCC)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(viridis)
library(ggh4x)
library(ggridges)

# Source code for plots
source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5

theme_set(theme_plot())

pal <- brewer.pal(name = "Dark2", n = 8)

Read and plot data

# Read from GH?
d <- read_csv("data/clean/trends.csv")
#> Rows: 92322 Columns: 17
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (3): species, life_stage, id
#> dbl (14): X, Y, quarter, temp_slope, biom_slope, phi_slope, mean_log_biomass...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(d)
#> Rows: 92,322
#> Columns: 17
#> $ X                   <dbl> 320.0167, 320.0167, 320.0167, 320.0167, 320.0167, …
#> $ Y                   <dbl> 6028.587, 6028.587, 6028.587, 6028.587, 6028.587, …
#> $ species             <chr> "cod", "cod", "cod", "cod", "dab", "dab", "dab", "…
#> $ life_stage          <chr> "adult", "adult", "juvenile", "juvenile", "adult",…
#> $ quarter             <dbl> 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1,…
#> $ temp_slope          <dbl> 0.03158410, 0.02985489, 0.03158410, 0.02985489, 0.…
#> $ biom_slope          <dbl> -0.5817831, -0.7105018, -0.4165413, -0.4676043, 4.…
#> $ phi_slope           <dbl> 0.001582802, 0.004173800, 0.004853177, 0.012797675…
#> $ id                  <chr> "320.016716551297_6028.58668110674_cod_adult_1", "…
#> $ mean_log_biomass    <dbl> 4.3916628, 4.5915365, 2.8446881, 2.9603249, 3.6929…
#> $ mean_temp           <dbl> 3.369871, 10.260329, 3.369871, 10.260329, 3.369871…
#> $ mean_oxy            <dbl> 8.114926, 6.209476, 8.114926, 6.209476, 8.114926, …
#> $ mean_phi            <dbl> 1.5023498, 0.8518787, 4.6064940, 2.6120242, 4.4595…
#> $ mean_log_biomass_ct <dbl> 1.544161015, 1.744034754, -0.002813649, 0.11282313…
#> $ mean_temp_ct        <dbl> -2.929083, 3.961374, -2.929083, 3.961374, -2.92908…
#> $ mean_oxy_ct         <dbl> 2.3680513, 0.4626013, 2.3680513, 0.4626013, 2.3680…
#> $ mean_phi_ct         <dbl> -1.6575600, -2.3080311, 1.4465842, -0.5478856, 1.2…

plot_map + 
  geom_raster(data = d, aes(X*1000, Y*1000, fill = biom_slope)) +
  facet_grid(life_stage ~ species) + 
  #scale_fill_gradient2() +
  scale_fill_viridis()
#> Warning: Removed 2488 rows containing missing values (geom_raster).


plot_map + 
  geom_raster(data = d, aes(X*1000, Y*1000, fill = temp_slope)) +
  facet_grid(life_stage ~ species) + 
  scale_fill_viridis()
#> Warning: Removed 2488 rows containing missing values (geom_raster).


plot_map + 
  geom_raster(data = d, aes(X*1000, Y*1000, fill = phi_slope)) +
  facet_grid(life_stage ~ species) + 
  scale_fill_viridis()
#> Warning: Removed 2488 rows containing missing values (geom_raster).


d_cod <- filter(d, species == "cod" & life_stage == "adult")
#> filter: removed 73,974 rows (80%), 18,348 rows remaining

p1 <- plot_map + 
  geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = biom_slope)) + 
  scale_fill_gradient2(name = "Biotic trend") + 
  geom_sf(size = 0.3)

p2 <- plot_map + 
  geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = temp_slope)) + 
  scale_fill_gradient2(name = "Temperature trend", breaks = c(0.025, 0.075, 0.125)) + 
  geom_sf(size = 0.3)

p3 <- plot_map + 
  geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = phi_slope)) + 
  scale_fill_gradient2(name = "\u03C6 trend") + 
  geom_sf(size = 0.3)

((p2 + p3) / (p1 + plot_spacer()))
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.

ggsave("figures/map_trends_cod.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted. Consider using geom_tile() instead.
#> Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.

# Scale biomass trends (temperature, phi)
d <- d %>% 
  mutate(group = paste(species, life_stage, sep = "_")) %>% 
  group_by(group) %>% 
  mutate(temp_slope_sc = (temp_slope - mean(temp_slope)) / sd(temp_slope),
         phi_slope_sc = (phi_slope - mean(phi_slope)) / sd(phi_slope)) %>% 
  ungroup() %>% 
  mutate(quarter_f = as.factor(quarter))
#> mutate: new variable 'group' (character) with 8 unique values and 0% NA
#> group_by: one grouping variable (group)
#> mutate (grouped): new variable 'temp_slope_sc' (double) with 92,310 unique values and 0% NA
#>                   new variable 'phi_slope_sc' (double) with 92,315 unique values and 0% NA
#> ungroup: no grouping variables
#> mutate: new variable 'quarter_f' (factor) with 2 unique values and 0% NA

Plot phi trend models

# Plot coefficients
phi_coefs %>% 
  filter(term %in% c("phi_slope_sc", "phi_slope_sc:mean_phi_ct")) %>% 
  ggplot(aes(estimate, group)) + 
  geom_point() + 
  geom_vline(xintercept = 0, linetype = 2, color = "grey30") + 
  facet_wrap(~term, ncol = 1)
#> filter: removed 32 rows (67%), 16 rows remaining


phi_pred_delta2 <- phi_pred_delta %>%
  # pivot_wider(names_from = mean_clim, values_from = bio_trend_delta) %>%
  # mutate(diff = high_phi - low_phi) %>%
  # mutate(Expected = ifelse(diff < 0, "No", "Yes")) %>%
  # We expect in general the biotic trends to be negative
  # We also expect biotic trends to be MORE negative if it's already in a low phi area
  #pivot_longer(c(low_phi, high_phi), names_to = "mean_clim", values_to = "bio_trend_delta") %>%
  mutate(mean_clim = ifelse(mean_clim == "high_phi", "High \u03C6", "Low \u03C6")) %>%
  mutate(id2 = paste(species, life_stage))
#> mutate: changed 16,000 values (100%) of 'mean_clim' (0 new NA)
#> mutate: new variable 'id2' (character) with 8 unique values and 0% NA
  
order <- phi_pred_delta2 %>%
  group_by(id2, mean_clim) %>% 
  summarise(mean_delta = mean(bio_trend_delta)) %>% 
  pivot_wider(names_from = mean_clim, values_from = mean_delta) %>% 
  mutate(diff = `Low φ` - `High φ`) %>% 
  arrange(desc(diff)) %>% 
  mutate(Expected = ifelse(`Low φ` < `High φ`, "Yes", "No"))
#> group_by: 2 grouping variables (id2, mean_clim)
#> summarise: now 16 rows and 3 columns, one group variable remaining (id2)
#> pivot_wider: reorganized (mean_clim, mean_delta) into (High φ, Low φ) [was 16x3, now 8x3]
#> mutate (grouped): new variable 'diff' (double) with 8 unique values and 0% NA
#> mutate (grouped): new variable 'Expected' (character) with 2 unique values and 0% NA

order
#> # A tibble: 8 × 5
#> # Groups:   id2 [8]
#>   id2               `High φ`  `Low φ`     diff Expected
#>   <chr>                <dbl>    <dbl>    <dbl> <chr>   
#> 1 Cod Adult          -1.08   -0.608    0.471   No      
#> 2 Cod Juvenile       -0.248   0.0758   0.323   No      
#> 3 Plaice Juvenile    -0.0106  0.00798  0.0185  No      
#> 4 Flounder Juvenile   0.0106  0.00898 -0.00158 Yes     
#> 5 Dab Juvenile        0.0366 -0.0206  -0.0572  Yes     
#> 6 Flounder Adult      0.233   0.0883  -0.144   Yes     
#> 7 Plaice Adult        0.0617 -0.204   -0.266   Yes     
#> 8 Dab Adult           0.404  -0.0714  -0.475   Yes

phi_pred_delta2 <- phi_pred_delta2 %>% 
  left_join(order %>% dplyr::select(Expected), by = "id2")
#> Adding missing grouping variables: `id2`
#> left_join: added one column (Expected)
#>            > rows only in x        0
#>            > rows only in y  (     0)
#>            > matched rows     16,000
#>            >                 ========
#>            > rows total       16,000

# phi_pred_delta2 %>%
#   ggplot(aes(bio_trend_delta, fct_reorder(id2, diff), color = mean_clim, alpha = Expected)) +
#   geom_vline(xintercept = 0, linetype = 2, color = "grey30") + 
#   geom_point(size = 3, position = position_dodge(width = 0.1)) +
#   scale_color_manual(values = c("steelblue2", "tomato3")) +
#   scale_alpha_manual(values = c(0.4, 1)) +
#   labs(y = "", x = "\u0394 biotic trend with a 1 st.dev\ndecline in \u03C6 trend", color = "Mean climate") +
#   theme_plot(base_size = 16) +
#   theme(legend.direction = "vertical") +
#   NULL

phi_pred_delta2$id2 <- factor(phi_pred_delta2$id2, levels = c(order$id2))
                                 
phi_pred_delta2 %>%
  ggplot(aes(bio_trend_delta, id2, fill = mean_clim, alpha = Expected)) +
  geom_vline(xintercept = 0, linetype = 2, color = "grey30") + 
  geom_density_ridges(color = NA) +
  scale_fill_manual(values = c("steelblue2", "tomato3")) +
  scale_alpha_manual(values = c(0.4, 0.8)) +
  labs(y = "", x = "\u0394 biotic trend with a 1 st.dev\ndecline in \u03C6 trend", fill = "Mean climate") +
  theme_plot(base_size = 16) +
  theme(legend.direction = "vertical") +
  NULL
#> Picking joint bandwidth of 0.00458


ggsave("figures/mi_trend_delta.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)
#> Picking joint bandwidth of 0.00458